
# Sanghyun Hong, February 23, 2019

#rm(list = ls())
###################################
## Simulation Environment Setup  ##
#################################################################
#pathname<-"C:/Temp/SimulationProject/"
pathnameAK<-paste(pathname,"SimCodes/AK_Codes/",sep="")
setwd(pathnameAK)
source("Modified_AllApplications.R")
alpha_List<-c(0.0,0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0);
Bias_List<-c('Sig');
Sim_Environ<-c('PRE')
SimN<-simN;
#################################################################
#
#
#
#
#
#################################################################
# Import Packages
#################################################################
library(sqldf) # This pakcage is for SQL
library(gridExtra)
library(grid)
library(lmtest)
library(sandwich)
library(metafor)
library(clubSandwich)
#################################################################
#
#
#
#
#
###############################################
## Primary Study Data (Effect) Generation    ##
#################################################################
PrimaryStudy <- function(StudyID, al, ali, lambdai0, type){
  if(type=='PRE'){
    sigr<-sqrt(0.25);
    m<-10;
    lambdai<-0.5+30*lambdai0;
  }else if(type=='RE'){
    sigr<-0;
    m<-1;
    lambdai<-0.5+30*lambdai0;
  }else if(type=='FE'){
    sigr<-0;
    m<-1;
    lambdai<-0.2+30*lambdai0;
  }
  obs<-100;
  PrimaryData<-as.data.frame(matrix(ncol=18, nrow=m))
  colnames(PrimaryData)<-c('StdID','EstID','al','ali','alir','lambdai','lambdair','effect','se','lowerBound','upperBound','Sig','PosSig','NegSig','nSig','nPosSig','nNegSig','nPos')
  for(i in 1:m){
    alir<-rnorm(1, mean = ali, sd = sigr);
    lambdair<- lambdai + runif(1, 0, 1)*(sigr>0)
    x <- rnorm(obs, mean = 0, sd = 1)
    y <- 1 + alir*x + lambdair*rnorm(obs, mean = 0, sd = 1)
    eff_se_pval<-as.numeric(summary(lm(y~x))$coefficients[2,c(1,2,4)])
    t<-abs(as.numeric(summary(lm(y~x))$coefficients[2,3]))
    ci<-as.numeric(confint(lm(y~x), level=0.95)[2,1:2])
    PrimaryData[i,1:14]<-c(StudyID, i, al,ali, alir, lambdai, lambdair, eff_se_pval[1:2], ci, (t>=2), (ci[1]>0), (ci[2]<0))
  }
  PrimaryData$nSig<-sum(PrimaryData$Sig)/m;
  PrimaryData$nPosSig<-sum(PrimaryData$PosSig)/m;
  PrimaryData$nNegSig<-sum(PrimaryData$NegSig)/m;
  PrimaryData$nPos<-sum(PrimaryData$effect>0)/m;
  return(PrimaryData)
}
###################################
## Meta Analysis Data Generation ##
#################################################################
MetaStudy <- function(type, alpha){
  if(type=='PRE'){
    StudyN<-100;
    sigi<-2
  }else if(type=='RE'){
    StudyN<-1000;
    sigi<-1
  }else if(type=='FE'){
    StudyN<-1000;
    sigi<-0;
  }
  
  for(i in 1:StudyN){
    ali<-rnorm(1, mean = alpha, sd = sigi)
    if(i==1){
      MetaStudyData<-PrimaryStudy(i, alpha, ali, runif(1,0,1), type)
    }else{
      MetaStudyData<-rbind(MetaStudyData, PrimaryStudy(i, alpha, ali, runif(1,0,1), type))
    }
  }
  return(MetaStudyData)
}
#################################################################
#################################################################
#
#
# 
#
#
#################################################################
# Clustered Standard Errors
#################################################################
clusteredSE <-function(regOLS, Study){
  M <- length(unique(Study)) 
  N <- length(Study) 
  K <- regOLS$rank
  dfc <- (M/(M-1)) * ((N-1)/(N-K))
  u<-apply(estfun(regOLS),2,function(x) tapply(x, Study,sum))
  vcovCL<-dfc*sandwich(regOLS, meat=crossprod(u)/N)
  ci <-coef(regOLS) + sqrt(diag(vcovCL)) %o% qt(c(0.025,0.975),summary(regOLS)$df[2])
  return (list("co"=coeftest(regOLS, vcovCL), "ci"=ci))
}
#################################################################
#
#
# 
#
#
#################################################
## Main Simulation For Table VII Begins Here   ##
#################################################################
start.time <- Sys.time();
nnnn<-length(Sim_Environ)*length(Bias_List)*length(alpha_List)*SimN
TableVII_Raw<-as.data.frame(matrix(NA, nrow=nnnn, ncol=25))
SubTable<-as.data.frame(matrix(NA, nrow=nnnn, ncol=16))
colnames(TableVII_Raw)<-c('Bias','Type','alpha','FPP_Bhat','FPP_Bias','FPP_Cov','WLSFE_Bhat','WLSFE_Bias','WLSFE_Cov','WLSRE_Bhat','WLSRE_Bias','WLSRE_Cov','WAAP_Bhat','WAAP_Bias','WAAP_Cov','WAAP2_Bhat','WAAP2_Bias','WAAP2_Cov','AK1_Bhat','AK1_Bias','AK1_Cov','AK3_Bhat','AK3_Bias','AK3_Cov',"I2")
colnames(SubTable)<-c('Bias','Type','alpha','FppCnt','WappCnt','WappTotal','WappSample','WappFrac','Wapp2Cnt','Wapp2Total','Wapp2Sample','Wapp2Frac','AK1P','AK3P1','AK3P2','AK3P3')
cnt<-0;
for(type in Sim_Environ){
  for(bias in Bias_List){
    for(alpha in alpha_List){
      for(i in 1:SimN){
        cnt<-cnt+1;
        MetaData<-as.data.frame(MetaStudy(type, alpha))
        MetaData<-cbind(MetaData, se2=MetaData$se*MetaData$se, invSE=1/MetaData$se, t=MetaData$effect/MetaData$se)
        MetaData<-as.data.frame(cbind(MetaData, rnd=runif(nrow(MetaData),0,1)))
          if(type=='PRE'){
          if(bias=='Sig'){
            MetaData<-subset(MetaData, (MetaData$nSig>=0.7));
          }else{
            MetaData<-subset(MetaData, (MetaData$nPos>=0.7));
          }
        }else{
          if(bias=='Sig'){
            MetaData<-subset(MetaData, ((MetaData$Sig==1)   | (MetaData$rnd<0.1)));
          }else{
            MetaData<-subset(MetaData, ((MetaData$effect>0) | (MetaData$rnd<0.1)));
          }
        }
        if(type=='PRE'){
          ###################
          ## PET-PEESE     ##
          ## PET-PEESE     ##
          #######################################################
          regTmp <- lm(t~1+invSE ,data=MetaData)
          regFPP <- clusteredSE(regTmp, MetaData$StdID)
          #coef_test(regTmp, vcov = "CR1", cluster = MetaData$StdID)
          fppcnt<-0
          if(regFPP$co[2,3]>1.96){
            regTmp <- lm(t~se+invSE-1 ,data=MetaData)
            regFPP <- clusteredSE(regTmp, MetaData$StdID)
            fppcnt<-1
          }
          FPP<-as.numeric(c(regFPP$co[2,1],regFPP$co[2,1]-alpha,(alpha>as.numeric(regFPP$ci[2,1]))*(alpha<as.numeric(regFPP$ci[2,2]))))
          #######################################################

          ###################
          ## WLS - FE      ##
          ## WLS - FE      ##
          #######################################################
          regTmp <- lm(t~invSE-1 ,data=MetaData)
          regFE <- clusteredSE(regTmp, MetaData$StdID)
          WLSFE<-as.numeric(c(regFE$co[1],regFE$co[1]-alpha,(alpha>as.numeric(regFE$ci[1]))*(alpha<as.numeric(regFE$ci[2]))))
          #######################################################

          ###################
          ## WLS - RE      ##
          ## WLS - RE      ##
          #######################################################
          rmaRE<-rma(effect~1, data=MetaData, vi=se*se, method="DL")
          tau2<-rmaRE$tau2;
          I2<-rmaRE$I2/100;
          re_weight<-1/(MetaData$se^2 + tau2)
          regTmp <- lm(effect~1, weights=re_weight, data=MetaData)
          regRE <- clusteredSE(regTmp, MetaData$StdID)
          WLSRE<-as.numeric(c(regRE$co[1],regRE$co[1]-alpha,(alpha>as.numeric(regRE$ci[1]))*(alpha<as.numeric(regRE$ci[2]))))
          #######################################################

          ###################
          ## WAAP          ##
          ## WAAP          ##
          #######################################################
          wappcnt<-1;
          effectWLS<-as.matrix(MetaData$effect)
          wWLS<-1/MetaData$se^2
          u_WLS<-sum(effectWLS*wWLS)/sum(wWLS)
          
          APSdata<-subset(MetaData, (abs(u_WLS/2.8)>=MetaData$se))
          if(nrow(APSdata)<2){APSdata<-MetaData; wappcnt<-0;}
          APSdata<-cbind(APSdata, u_WLS=u_WLS)
          
          wapptot <- nrow(MetaData);
          wappsam <- nrow(APSdata);
          wappfrac <- wappsam/wapptot;
          
          regTmp <- lm(t~invSE-1 ,data=APSdata)
          regWAAP <- clusteredSE(regTmp, APSdata$StdID)
          WAAP<-as.numeric(c(regWAAP$co[1],regWAAP$co[1]-alpha,(alpha>as.numeric(regWAAP$ci[1]))*(alpha<as.numeric(regWAAP$ci[2]))))
          #######################################################

          ###################
          ## WAAP2         ##
          ## WAAP2         ##
          #######################################################
          wapp2cnt<-1;
          tstat_waap2<-MetaData$effect/MetaData$se;
          invSE_waap2<-1/MetaData$se;
          u_WLS<-as.numeric(coefficients(lm(tstat_waap2~invSE_waap2))[2]);
          
          APSdata<-subset(MetaData, (abs(u_WLS/2.8)>=MetaData$se))
          if(nrow(APSdata)<2){APSdata<-MetaData; wapp2cnt<-0;}
          APSdata<-cbind(APSdata, u_WLS=u_WLS)
          
          wapp2tot <- nrow(MetaData);
          wapp2sam <- nrow(APSdata);
          wapp2frac <- wapp2sam/wapp2tot;
          regTmp <- lm(t~invSE-1 ,data=APSdata)
          regWAAP <- clusteredSE(regTmp, APSdata$StdID)
          WAAP2<-as.numeric(c(regWAAP$co[1],regWAAP$co[1]-alpha,(alpha>as.numeric(regWAAP$ci[1]))*(alpha<as.numeric(regWAAP$ci[2]))))
          #######################################################
        }
        if(type=='RE'){
          ###################
          ## PET-PEESE     ##
          ## PET-PEESE     ##
          #######################################################
          fppcnt<-0
          regTmp <- lm(t~1+invSE ,data=MetaData)
          regFPP <- coeftest(regTmp, vcov = vcovHC(regTmp, "HC1"))[1:2, c(1,3,4)]
          regci <- confint(regTmp, vcov = vcovHC(regTmp, "HC1"))[2,1:2]
          
          if(regFPP[2,2]>1.96){
            fppcnt<-1
            regTmp <- lm(t~se+invSE-1 ,data=MetaData)
            regFPP <- coeftest(regTmp, vcov = vcovHC(regTmp, "HC1"))[1:2, c(1,3,4)]
            regci <- confint(regTmp, vcov = vcovHC(regTmp, "HC1"))[2,1:2]
          }
          FPP<-as.numeric(c(regFPP[2,1],regFPP[2,1]-alpha,(alpha>as.numeric(regci[1]))*(alpha<as.numeric(regci[2]))))
          #######################################################

          ###################
          ## WLS - FE      ##
          ## WLS - FE      ##
          #######################################################
          regTmp <- lm(t~invSE-1 ,data=MetaData)
          regFE <- as.numeric(coeftest(regTmp, vcov = vcovHC(regTmp, "HC1"))[1])
          regci <- as.numeric(confint(regTmp, vcov = vcovHC(regTmp, "HC1"))[1:2])
          WLSFE<-as.numeric(c(regFE,regFE-alpha,(alpha>regci[1])*(alpha<regci[2])))
          #######################################################

          ###################
          ## WLS - RE      ##
          ## WLS - RE      ##
          #######################################################
          rmaRE<-rma(effect~1, data=MetaData, vi=se*se, method="DL")
          tau2<-rmaRE$tau2;
          I2<-rmaRE$I2/100;
          re_weight<-1/(MetaData$se^2 + tau2)
          regTmp <- lm(effect~1,data=MetaData, weights = re_weight)
          regRE <- as.numeric(coeftest(regTmp, vcov = vcovHC(regTmp, "HC1"))[1])
          regci <- as.numeric(confint(regTmp, vcov = vcovHC(regTmp, "HC1"))[1:2])
          WLSRE<-as.numeric(c(regRE,regRE-alpha,(alpha>regci[1])*(alpha<regci[2])))
          #######################################################

          ###################
          ## WAAP          ##
          ## WAAP          ##
          #######################################################
          wappcnt<-1;
          effectWLS<-as.matrix(MetaData$effect)
          wWLS<-1/MetaData$se^2
          u_WLS<-sum(effectWLS*wWLS)/sum(wWLS)
          APSdata<-subset(MetaData, (abs(u_WLS/2.8)>=MetaData$se))
          if(nrow(APSdata)<2){APSdata<-MetaData; wappcnt<-0;}
          wapptot <- nrow(MetaData);
          wappsam <- nrow(APSdata);
          wappfrac <- wappsam/wapptot;
          
          regTmp <- lm(t~invSE-1 ,data=APSdata)
          regWAAP <- as.numeric(coeftest(regTmp, vcov = vcovHC(regTmp, "HC1"))[1])
          regci <- as.numeric(confint(regTmp, vcov = vcovHC(regTmp, "HC1"))[1:2])
          WAAP<-as.numeric(c(regWAAP,regWAAP-alpha,(alpha>regci[1])*(alpha<regci[2])))
          #######################################################

          ###################
          ## WAAP2         ##
          ## WAAP2         ##
          #######################################################
          wapp2cnt<-1;
          tstat_waap2<-MetaData$effect/MetaData$se;
          invSE_waap2<-1/MetaData$se;
          u_WLS<-as.numeric(coefficients(lm(tstat_waap2~invSE_waap2))[2]);
          
          APSdata<-subset(MetaData, (abs(u_WLS/2.8)>=MetaData$se))
          if(nrow(APSdata)<2){APSdata<-MetaData; wapp2cnt<-0;}
          APSdata<-cbind(APSdata, u_WLS=u_WLS)
          
          wapp2tot <- nrow(MetaData);
          wapp2sam <- nrow(APSdata);
          wapp2frac <- wapp2sam/wapp2tot;
          
          regTmp <- lm(t~invSE-1 ,data=APSdata)
          regWAAP2 <- as.numeric(coeftest(regTmp, vcov = vcovHC(regTmp, "HC1"))[1])
          regci <- as.numeric(confint(regTmp, vcov = vcovHC(regTmp, "HC1"))[1:2])
          WAAP2<-as.numeric(c(regWAAP2,regWAAP2-alpha,(alpha>regci[1])*(alpha<regci[2])))
          #######################################################
        }
          
        ###################
        ## A-K           ##
        ## A-K           ##
        #######################################################
        AKdata<-as.data.frame(cbind(MetaData$effect, MetaData$se, MetaData$StdID))
        attach(setup(AKdata, TRUE), warn.conflicts = FALSE)
        source("EstimatingSelection.R")
        ciAK<-as.numeric(Psihat[1]+qt(c(.025, .975), df=(nrow(my.data)-1))*se_robust[1])
        AK1<-c(Psihat[1], Psihat[1]-alpha, (alpha>ciAK[1])*(alpha<ciAK[2]))
        Ak1p1<-Psihat[3]
        detach("setup(AKdata, TRUE)")
        #######################################################

        ###################
        ## A-K 3         ##
        ## A-K 3         ##
        #######################################################
        AKdata<-as.data.frame(cbind(MetaData$effect, MetaData$se, MetaData$StdID))
        attach(setup(AKdata, FALSE), warn.conflicts = FALSE)
        source("EstimatingSelection.R")
        ciAK<-as.numeric(Psihat[1]+qt(c(.025, .975), df=(nrow(my.data)-1))*se_robust[1])
        AK3<-c(Psihat[1], Psihat[1]-alpha, (alpha>ciAK[1])*(alpha<ciAK[2]))
        Ak3p1<-Psihat[3]; Ak3p2<-Psihat[4]; Ak3p3<-Psihat[5];
        detach("setup(AKdata, FALSE)")
        #AK3<-c(0,0,0)
        #######################################################
        TableVII_Raw[cnt,]<-c(bias,type,alpha, FPP, WLSFE, WLSRE, WAAP, WAAP2, AK1, AK3, I2)
        SubTable[cnt,]<-c(bias,type,alpha, fppcnt, wappcnt, wapptot, wappsam, wappfrac, wapp2cnt, wapp2tot, wapp2sam, wapp2frac, Ak1p1, Ak3p1, Ak3p2, Ak3p3)
                          
        cat("Simulation Environment #6: ", cnt, "/", nnnn,"\n");
      }
}}}
TableVII_Raw$FPP_Bhat<-as.numeric(TableVII_Raw$FPP_Bhat)
TableVII_Raw$FPP_Bias<-as.numeric(TableVII_Raw$FPP_Bias)
TableVII_Raw$FPP_Cov<-as.numeric(TableVII_Raw$FPP_Cov)

TableVII_Raw$WLSFE_Bhat<-as.numeric(TableVII_Raw$WLSFE_Bhat)
TableVII_Raw$WLSFE_Bias<-as.numeric(TableVII_Raw$WLSFE_Bias)
TableVII_Raw$WLSFE_Cov<-as.numeric(TableVII_Raw$WLSFE_Cov)

TableVII_Raw$WLSRE_Bhat<-as.numeric(TableVII_Raw$WLSRE_Bhat)
TableVII_Raw$WLSRE_Bias<-as.numeric(TableVII_Raw$WLSRE_Bias)
TableVII_Raw$WLSRE_Cov<-as.numeric(TableVII_Raw$WLSRE_Cov)

TableVII_Raw$WAAP_Bhat<-as.numeric(TableVII_Raw$WAAP_Bhat)
TableVII_Raw$WAAP_Bias<-as.numeric(TableVII_Raw$WAAP_Bias)
TableVII_Raw$WAAP_Cov<-as.numeric(TableVII_Raw$WAAP_Cov)

TableVII_Raw$WAAP2_Bhat<-as.numeric(TableVII_Raw$WAAP2_Bhat)
TableVII_Raw$WAAP2_Bias<-as.numeric(TableVII_Raw$WAAP2_Bias)
TableVII_Raw$WAAP2_Cov<-as.numeric(TableVII_Raw$WAAP2_Cov)

TableVII_Raw$AK1_Bhat<-as.numeric(TableVII_Raw$AK1_Bhat)
TableVII_Raw$AK1_Bias<-as.numeric(TableVII_Raw$AK1_Bias)
TableVII_Raw$AK1_Cov<-as.numeric(TableVII_Raw$AK1_Cov)

TableVII_Raw$AK3_Bhat<-as.numeric(TableVII_Raw$AK3_Bhat)
TableVII_Raw$AK3_Bias<-as.numeric(TableVII_Raw$AK3_Bias)
TableVII_Raw$AK3_Cov<-as.numeric(TableVII_Raw$AK3_Cov)

TableVII_Raw$I2<-as.numeric(TableVII_Raw$I2)
TableVII<-sqldf('select Bias, Type, alpha, avg(FPP_Bhat) as FPP_Bhat, avg(FPP_Bias*FPP_Bias) as FPP_MSE, avg(FPP_Cov) as FPP_TypeI, avg(WLSFE_Bhat) as WLSFE_Bhat, avg(WLSFE_Bias*WLSFE_Bias) as WLSFE_MSE, avg(WLSFE_Cov) as WLSFE_TypeI, avg(WLSRE_Bhat) as WLSRE_Bhat, avg(WLSRE_Bias*WLSRE_Bias) as WLSRE_MSE, avg(WLSRE_Cov) as WLSRE_TypeI, avg(WAAP_Bhat) as WAAP_Bhat, avg(WAAP_Bias*WAAP_Bias) as WAAP_MSE, avg(WAAP_Cov) as WAAP_TypeI, avg(WAAP2_Bhat) as WAAP2_Bhat, avg(WAAP2_Bias*WAAP2_Bias) as WAAP2_MSE, avg(WAAP2_Cov) as WAAP2_TypeI, avg(AK1_Bhat) as AK1_Bhat, avg(AK1_Bias*AK1_Bias) as AK1_MSE, avg(AK1_Cov) as AK1_TypeI , avg(AK3_Bhat) as AK3_Bhat, avg(AK3_Bias*AK3_Bias) as AK3_MSE, avg(AK3_Cov) as AK3_TypeI, avg(I2) as I2  from TableVII_Raw group by Type,Bias,alpha')
SubTable_Ave<-sqldf('select Bias, Type, alpha, avg(FppCnt) as FppCnt, avg(WappCnt) as WappCnt, avg(WappTotal) as WappTotal, avg(WappSample) as WappSample, avg(WappFrac) as WappFrac, avg(Wapp2Cnt) as Wapp2Cnt, avg(Wapp2Total) as Wapp2Total, avg(Wapp2Sample) as Wapp2Sample, avg(Wapp2Frac) as Wapp2Frac, avg(AK1P) as AK1P, avg(AK3P1) as AK3P1, avg(AK3P2) as AK3P2, avg(AK3P3) as AK3P3 from SubTable group by Type,Bias,alpha')

write.csv(TableVII, paste(pathname,"/Output/SimulationEnvironment06.csv", sep=""))
write.csv(SubTable_Ave, paste(pathname,"/Output/SimulationEnvironment06_Miscellaneous.csv", sep=""))

end.time <- Sys.time();
time.taken <- end.time - start.time;
cat("Total Time Taken: ", time.taken, "\n\n\n")
#################################################################
#################################################################